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The curvature dependence of the surface tension is related to the excess equimolar radius of 
liquid drops, i.e., the deviation of the equimolar radius from that defined with the macroscopic 
capillarity approximation. Based on the Tolman [J. Chem. Phys. 17, 333 (1949)] approach and its 
interpretation by Nijmeijer et al. [J. Chem. Phys. 96, 565 (1991)], the surface tension of spherical 
interfaces is analysed in terms of the pressure difference due to curvature. In the present study, the 
excess equimolar radius, which can be obtained directly from the density profile, is used instead of 
the Tolman length. Liquid drops of the truncated-shifted Lennard- Jones fiuid are investigated by 
molecular dynamics simulation in the canonical ensemble, with equimolar radii ranging from 4 to 
33 times the Lennard- Jones size parameter a. In these simulations, the magnitudes of the excess 
equimolar radius and the Tolman length are shown to be smaller than a/2. Other methodical 
approaches, from which mutually contradicting findings have been reported, are critically discussed, 
outlining possible sources of inaccuracy. 
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I. INTRODUCTION 

The macroscopic capillarity approximation consists in 
neglecting the curvature dependence of the surface ten- 
sion 7 of a spherical liquid drop. Accordingly, the surface 
teirsioir of a curved interface in equilibrium is approxim- 
ated by the value 70 in the zero-curvature limit, i.e., for a 
planar vapour-liquid interface. The Young-Laplace equa- 
tion pj43|] for spherical interfaces relates the macroscopic 
surface tension to a characteristic radius Rk of the liquid 
drop 

^ = lip'-p")^V, (1) 

which will be referred to as the capillarity radius here. 
Both the factor l/i?K and the difference between the li- 
quid pressure p' and the vapour pressure -p" characterize 
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the extent by which the surface is curved; the notation 
(/? = {p' —p")/2 for half of the pressure difference is intro- 
duced for convenience. At equilibrium, the temperature 
is the same for both phases, and the pressures p' and 
p" correspond to states with the same chemical poten- 
tial. The surface tension 70 of the planar vapour-liquid 
phase boundary, which is relatively easy to access exper- 
imentally, couples these two measures of curvature as a 
proportionality constant. 

In combination with an equation of state for the bulk 
fluid, microscopic properties such as the radius of a small 
liquid drop can thus be deduced from the macroscopic 
state of the surrounding vapour, i.e., from its supersat- 
uration ratio, and vice versa. This approach is the most 
widespread interpretation of the Gibbs theory of inter- 
faces [3, and it is the point of departure for the clas- 
sical nucleation theory (CNT) as introduced by Volmer 
and Weber [6] and further developed by Farkas [7j as 
well as subsequent authors [sl-fioj. With the Gibbs ap- 
proach one presumes a sharp dividing surface between 
the phases, a conceptual picture that does not reflect 
the physical phenomena present at the molecular length 
scale. However, this abstraction is precisely its strength. 
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Instead of discussing thermodynamic properties such as 
the density, the pressure tensor and the free energy dens- 
ity in a locaHzed way, interfacial excess quantities can be 
assigned to the formal dividing surface whole. 

It should be recognized that significant size effects on 
interfacial properties had already been detected experi- 
mentally by Weber [ll| at the turn of the last century. 
This was also known to Farkas 01 who stated explicitly 
that the capillarity approximation should be expected 
to fail for radii at the length scale of the intermolecular 
interactions. In the absence of a better approximation, 
however, the surface tension of the planar phase bound- 
ary had to be used for nucleation theory, and little has 
changed in this respect in the meantime. 

In case of significant deviations from the macroscopic 
capillarity approximation, liquid drops cannot be charac- 
terized sufficiently by a single effective radius. Instead, 
the capillarity radius i?^ is distinct from the equimolar 
radius Rp, which is also known as the Gibbs adsorp- 
tion radius. For a single-component system, the latter 
is defined by the zero excess density criterion 



dzz^ [p{z) - p'{p,T)] 

i 

dzz^ [p(z)-p"(Ai,r)]=0, 



(2) 



i.e., by comparing a step function based on the bulk li- 
quid and vapour number densities p'(/x,T) and p"(^,T) 
as functions of the chemical potential fi and the temper- 
ature T, respectively, with the microscopic radial density 
profile p{z). By convention, the density p corresponds 
to the number of particles per volume here, rather than 
their mass, and z denotes the distance from the centre 
of mass of the liquid drop. In the following discussion, 
T is treated as a parameter (instead of a variable), so 
that total differentials are to be understood as partial 
differentials at constant temperature. 

For curved interfaces in equilibrium, the chemical po- 
tential deviates from its saturated value /ig for a flat in- 
terface. In case of a drop, both phases are supersat- 
urated. To realize this, it is sufficient to consider the 
Gibbs-Duhem equation for a curved phase boundary 



d {p - v") = {p - p") dp.- 



(3) 



For a planar interface, both phases coexist at the satur- 
ation condition (/i = /ig) and the pressure difference is 
zero. Raising the value of the liquid pressure p' over the 
vapour pressure p" therefore increases the chemical po- 
tential fi, which must be equal for both phases in (stable 
or unstable) equilibrium, so that its value for a system 
with a liquid drop will exceed /ig- The precise conditions 
can be determined from the pressure difference between 
the fluid phases by means of an equation of state. 

Beside and Rp, a thermodynamically relevant defin- 
ition of the liquid drop size is given by the surface of 
tension radius 

1 



R'^ 



(4) 



which is also known as the Laplace radius. It can be ob- 
tained by inserting the actual value of the surface tension 
7 of the system with the curved interface (not the planar 
limit value) into the Young-Laplace equation. This ra- 
dius can be related to the surface area a and to the 
volume V of the drop 



R^da^2 dV. 



(5) 



The excess grand potential S of the surface thus evalu- 
ates to 



R^dS = 27 dy. 



in terms of the surface tension 



7 



dS 
da 



(6) 



(7) 



Modified versions of the Young-Laplace equation, which 
allow for the use of different radii in an analogous way, 
were introduced by Buff and Kondo [14.]. 

The present study deals with the deviation between 
the capillarity radius R^,, the equimolar radius Rp and 
the surface of tension radius R^ of a liquid drop in equi- 
librium with a supersaturated vapour. As Tolman [l5|- 
il7]] . following Gibbs, showed on the basis of axiomatic 
thermodynamics, one of these differences, now commonly 
referred to as the Tolman length 



Rp R^ , 



(8) 



is sufficient to characterize the curvature dependence of 
the surface tension Uji 
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It is important to point out that this relation is exact, 
strictly following the approach of Gibbs, i.e., without neg- 
lecting any of the higher-order curvature terms. The cu- 
bic expression derives from an integral over the spherical 
density profile. However, Eq. ^ is often transformed 
into a polynomial expansion for 70/7, which contains an 
infinite number of terms and has to be truncated, e.g., 
after the second-order contribution in terms of curvature 

m 



2^^ = 1 

7 



R 



(10) 



Here, 5o is the Tolman length in the zero-curvature (in- 
finite radius) limit. Castellanos et al. jl^ have conjec- 
tured that «the Tolman length is related to the interfacial 
width A'' according to A'^ w 25. » The Block length A, 
which characterizes the effect of Gaussian curvature that 
becomes predominant when 5 is very small or for systems 
where, due to an inherent symmetry, 5 = holds by con- 
struction, has recently been investigated by Block et al. 
[isj : a similar leading term, proportional to R^^lnR^, 
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has also been deduced by Bieker and Dietrich [2^ from 
DFT based on a Barker-Henderson perturbation expan- 
sion. 

One should keep in mind that the Tolman equation 
as given by Eq. ([9]) is valid for curved phase boundaries 
of pure fluids in general, whereas truncated polynomial 
expansions in terms of the curvature 1 / like Eq. (jlOp 
necessarily break down for liquid drops at the molecu- 
lar length scale. In practice, one of the major problems 
of the Tolman approach is that it analyses the surface 



tension in terms of the radii Rp and R^ . 



While Rp can 



be immediately obtained from the density profile, R^ is 
by definition related to 7 itself. Since for highly curved 
interfaces the value of 7 is disputed or unknown [2l|, [l^l , 
the surface of tension radius R-^ is correspondingly un- 
certain. 

To resolve this issue, we reformulate Tolman's theory 
in terms of R^ and Rp. This leads to greater transpar- 
ency, since the capillarity radius i?^ can be obtained on 
the basis of the surface tension in the planar limit 70, 
which is experimentally accessible, and properties of the 
(stable and metastable) bulk fluid. It is related to the 
pressure difference between the coexisting phases in equi- 
librium, which is a bulk property as well, since it can be 
determined from ^ and T with an equation of state for 
the fluid. All information on the molecular structure of 
the curved interface can thus be captured by a single 
undisputed quantity here, namely the equimolar radius 
Rp. 

For this approach, the excess equimolar radius, defined 



(11) 



plays a role similar to the Tolman length, and the mac- 
roscopic quantity ip is used instead of 1/ R^ as a measure 
of the influence of curvature on the thermophysical prop- 
erties of the interface and the bulk phases. In this way, 
the thermodynamics of liquid drops are discussed by fol- 
lowing a new route that relies on the density profiles and 
bulk properties only, avoiding the intricacies of defining 
the pressure tensor or the change in the surface area as 
required by other approaches. 

The present method is related to the «direct determin- 
ation» of Sq proposed by Nijmeijer et al. as recently 
applied by van Giessen and Blokhuis [24] on the basis of 
a representation of ipRp over 1/Rp with 
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lim — - — ■ 

Rp^oo d(l/Rp) 



^R, 



(12) 



as depicted in Fig. [T] However, the implementation sug- 
gested here is methodologically different from that of van 
Giessen and Blokhuis which relies on a pressure tensor to 
obtain if, whereas in the present work, the pressure dif- 
ference is determined by molecular dynamics (MD) sim- 
ulation of the bulk fiuids. Applying the definitions of the 
capillarity radius and the excess equimolar radius, Eq. 
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Figure 1. Representation of van Giessen and Blokhuis (23 |. 
showing ipRp as a function of the equimolar curvature 1/Rp 
for liquid drops of the truncated-shifted Lennard- Jones fluid 
at T = 0.9 e/k, where the equimolar radius Rp is determined 
from the density profiles and (/? from the difference between 
the values of the normal component of the Irving-Kirkwood 
pressure tensor in the homogeneous regions inside the liquid 
drop as well as outside, i.e., in the homogeneous supersatur- 
ated vapour. In comparison with the results of van Giessen 
and Blokhuis (□), the da,ta of Vrabec et al. 'i-")| (o), which 
were obtained by the same method, are included here along 
with a data point (A) where ip is determined by MD sim- 
ulation of the homogeneous fluid. The data for the planar 
surface tension 70 are taken from simulations of Vrabec et al. 
(•) and van Giessen and Blokhuis (■) as well as the correla- 
tion of Vrabec et al. (a). The continuous lines are guides to 
the eye: In the planar limit, a positive slope corresponds to a 
negative Tolman length and vice versa, cf. Eq. (|12p . 
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(13) 



facilitating an analysis of interface properties in terms of 
the radii R^ and Rp as well as the deviation -q between 
them. 

This article is structured as follows: In Section|lTl a re- 
view is made of the available routes to the Tolman length 
and the surface tension by molecular simulation. MD 
simulation methods immediately related to nucleation it- 
self, from which information of the excess free energy of 
curved interfaces can also be deduced f26l-[29l|. are not 
included in that discussion; in this regard, the reader is 
referred to Chkonia et al. [30|. Section [III] is dedicated to 
a brief outline of how Tolman's thermodynamic approach 
can be transformed by analysing the surface tension in 
terms of 77 and ip rather than 5 and 1/R^. The meth- 
odology and the results of a series of canonical ensemble 
MD simulations, where the excess equimolar radius is ob- 
tained solely on the basis of density profiles, are presented 
in Section IIVI An interpretation of these results is given 
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in Section |Vl placing the present findings in the context 
of the multitude of mutually contradicting hypotheses 
proposed in the literature. 



II. THE TOLMAN LENGTH FROM 
MOLECULAR SIMULATION 

A. Analysis of the planar interface 

For the planar interface, the definition of the Tol- 
man length given by Eq. (jS]) ceases to be applicable, as 
the surface of tension radius becomes ill-defined in 
the absence of curvature, since the pressure is equal on 
both sides of the interface in this case. Therefore, the 
planar interface Tolman length Sq necessarily has to be 
derived from considerations pertaining to curved geomet- 
ries. It can be obtained either by extrapolating results 
for d to the macroscopic limit iy9 — > (i.e., R-y — ^ oo) or 
by constructing the limit explicitly from expressions for 
the radii Rp and R^. The latter approach was followed 
by Fisher and Wortis [lH who, on the basis of Landau 
(square-gradient) theory, derived the relation 



mm 



So 



dpo{z) 



dpo{z] 



Z — — OC 



dz 
dpo{z) z, 



dpo(z) 



dpojz) 
dlnz 

(14) 



in terms of the density profile po{z) of the planar inter- 
face. This expression can also be extended to account 
for the pair density profile, whereby Eq. becomes a 
limiting case [H, [s^ . 

The available computational methods for evaluating 
the Tolman length of curved interfaces, however, involve 
the determination of the surface tension 7. It is usually 
the methodology related to the evaluation of 7 that is 
both the crucial and the most debatable step, which is 
made evident by the contradictory findings for 7 (and 
consequently also for Sq) obtained from different meth- 
ods. Three routes to the surface tension of liquid drops 
will now be discussed briefly: the mechanical route as 
implemented by Thompson et at. [s^l, the grand canon- 
ical route of Schrader et al. |35| and the variational route 
developed by Sampayo et al. |22l |. 

Many different versions and combinations of these ap- 
proaches exist (sgI - IssI . but it would be inappropriate to 
attempt a full appreciation of the complete body of work 
here. The reader is directed to the excellent review by 
Henderson [so'] for a detailed discussion of the underlying 
statistical mechanical approaches. 



B. The mechanical route 

The mechanical route to the surface tension is based 
on the Bakker-Buff equation for spherical interfaces (Tsl . 



7 



dzz"^ [Pn{z) -Pt{z)] 



(15) 



z=0 



in terms of the normal component Pn{z) and the two 
(equal) tangential components pt{z) of the diagonalized 
pressure tensor, which is considered as a spherical aver- 
age, and where the integration is from the centre of the 
drop (z= 0). With this relation one expresses the work 
required for a reversible isothermal deformation of the 
system that leads to an infinitesimal increase of the sur- 
face area at constant volume, which coincides with the 
associated free energy difference. It is sufficient to com- 
pute either the normal or the tangential pressure profile, 
since both are related by [13, 113 



dpn 
dlnz 



2 {Pt - Pn) 



(16) 



At mechanical equilibrium, Eq. (jl5|) can thus be trans- 
formed to [H 



27^ 



dpn{z) 



(17) 



z=0 



a term in which R^ no longer appears. The surface of ten- 
sion radius R-y can be obtained from the Young-Laplace 
equation once the surface tension 7 is known. 

The most widespread implementation of this approach 
in terms of intermolecular pair potentials makes use of 
the Irving-Kirkwood (IK) [43| pressure tensor, \vhich was 
first applied to (spherical) interfaces by Buff |l3| and un- 
derlies the simulation studies of Vrabec et al. [23 as well 
as those of van Giessen and Blokhuis [IJ]. Its normal 
component is given by [13, SI] 



p„(z) = fcrp(z)+ >' 



E 



drij Attz^ rij 



(18) 



wherein k is the Boltzmann constant and the summation 
covers the set S containing all sets of particles i and j that 
are connected by a line intersecting a sphere of radius z 
around the centre of mass. The intersection coordinates 
relative to the centre of mass of the liquid drop are rep- 
resented by z and the distance between the particles by 
Tij with rij = \rij\, while —duij/drij is the force acting 
between the two particles i and j. 

Regarding the mechanical route as described here, vari- 
ous issues arise: 

• It is not clear to what extent the spherical average 
of the pressure tensor succeeds in accounting for the 
free energy contribution of capillary waves, i.e., the 
excited vibrational modi of the interface (44l. |45|. 

• Irving and Kirkwood [43| originally proposed their 
expression for the special case of «a single compon- 
ent, single phase system». Its derivation relies on 
truncating an expansion in terms of derivatives of 
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the pair density p*^^) after the first term, thereby 
disregarding the density gradient completely. For 
a liquid drop, this can lead to inaccuracies: «at a 
boundary or interface . . . neglecting terms beyond 
the first may not be justified» [43|. 

• By construction, the mechanical route cannot be 
separated from the assumption of a mechanical 
equilibrium that underlies both the basic approach, 
i.e., Eqs. p5|) to (fT6| . and the derivation of the IK 
pressure tensor, cf. Eq. ([T^ . For nanoscopic hquid 
drops, however, configurations deviating from the 
equilibrium shape correspond to a significant frac- 
tion of the partition function. 

• The non-unique nature of the pressure tensor, 
which for a planar interface does not have a con- 
sequence on the computed value of the surface 
tension [i^, leads to an inconsistent description 
for a curved interface (ssl . Isol . l47l |. However, the 
Harasima pressure tensor [41], where the set S is 
defined differently and the tangential pressure pro- 
file pt{z) is computed instead of the normal com- 
ponent Pn{z), has been found to agree rather well 
with the IK tensor f2l,[33,|4i. 



C. The grand canonical route 

Prom an analysis of the canonical partition function 
and its dependence on the characteristic length L of 
otherwise similar systems. Binder (45j derived very use- 
ful scaling laws for the probability Lu{pmin) of relatively 
small subvolume to have the density pmin corresponding 
to a maximum of the local free energy, i.e., the least prob- 
able local density between p' and p" . It follows that «the 
probability of a homogeneous state with order parameter 
Pmin decreases exponentially fast with the volume» while 
for cases where the corresponding subvolume is situated 
within a phase boundary the probabilit y << decreases ex- 
ponentially fast with the interface area» [43| . The surface 
excess of the grand potential (per unit surface area) can 
thus be determined as 



— lim — — kT lim ttT" 

a— s-oo a L—KX o,\L) 



(19) 



which is related to the surface tension by 7 = dS /da. 
Therein, the term a{L) describes the dependence of the 
surface area on the characteristic length of the system 
[isj . e.g., a{L) = 2L^ for a planar slab in a cubic volume 
V = with standard periodic boundary conditions. 

Small subvolumes of a canonical system in the thermo- 
dynamic limit (iV —^00) are equivalent to systems with 
constant /x, V and T so that grand canonical Monte Carlo 
(GCMC) simulation can equally be applied. Umbrella 
sampling may be used to fully sample the relevant range 
of values for the order parameter [3 E^l , corresponding 
to the number of particles N present in the grand canon- 
ical system. Thereby, a profile is obtained for the free 



energy density f[N) or /(p), i.e., the dependence of the 
grand potential per volume unit on the order parameter. 

To analyse liquid drops of a certain size, however, the 
limit a — > 00 cannot be applied since the area a of the 
surface of tension is fixed. Instead, the surface excess 
term f^(Rp) is determined from expressions based on 
the equimolar radius (35l | 



fip) = yp'iRpMRp) 
+ ^p"{Rp)p{Rp) 



(20) 



V' = AnRp/S is the volume associated with the liquid 
phase here, V" — V — V is the remainder of the volume 
and p'{Rp) as well as p"{Rp) are bulk densities related to 
the liquid drop and the surrounding vapour. The chem- 
ical potential p{Rp) is equal for the vapour and liquid 
regions, but different from both the saturated bulk value 
/is and the chemical potential p used for the grand canon- 
ical simulation itself. This formalism has recently been 
employed by Schrader et al. [2ll . [35| as well as Block et 
al. [18], to examine the interfacial properties of drops, 
bubbles, and symmetric mixtures in great detail. 

The original method of Binder [45] was developed for 
planar interfaces. In the case of systems with a spherical 
geometry, the following points should be kept in mind: 

• Following the approach of Schrader et al. [sH] , the 
surface tension 7 can be accessed only indirectly, 
e.g., from Eq. based on the surface of tension 
radius R-y which also has to be obtained in a cir- 
cuitous manner. Thereby, care should be taken not 
to confuse with 7, or Rp with R^. 

• Since the infinite size limit, cf. Eq. ([T^ . does not 
apply to nanoscopic liquid drops and the systems 
under consideration can be extremely small, it is 
not generally possible to neglect the contribution 
of homogeneous configurations to f{p) [45| . 

• The assumption that AttR^ is the surface area as- 
sociated with the surface excess for the grand po- 
tential of the system, as in Eq. (PU)) . essentially 
amounts to applying the macroscopic capillarity 
approximation. Such an approach may be justified 
under certain circumstances, but for investigations 
of the deviation from capillarity it is of limited use 
only. 

Other umbrella sampling based methods (sol. [5l|. which 
will not be discussed here in detail, are confronted with 
similar difficulties, in particular regarding the relation 
between the surface tension and the surface excess free 
energy. 



D. The variational route 

The variational route to the surface tension is based 
on Bennett's [52,] general considerations of the molecular 



6 



simulation of free energies and entropic quantites. In 
the canonical ensemble, the free energy difference A A — 
Ai — Aq between two states with equal N, V and T is 
given by the quotient of the respective canonical partition 
functions Zq and Zi , which can be evaluated as averages 
in terms of internal energy differences (52| 



exp 



{AA\ ^ Zo 



V kT 



Zi 

(min (1, exp {[Ei 



E,]/[kT]))), 



,(21) 



where the index of the angular brackets denotes the sys- 
tem over which an ensemble average is taken. Bennett 
proposed the determination of these energy differences 
from «separately-generated samples» [s^l for Eq and Ei . 
If the two systems differ in the area of a phase boundary, 
then the free energy difference can be related to the sur- 
face tension, assuming that all other deviations between 
the two states are accurately taken into account. 

Gloor et al. [3fl| introduced a version of this approach 
where differences between the two states are obtained 
from a single simulation run for an unperturbed system 
with the partition function Zq. Corresponding config- 
urations of the second, perturbed system are generated 
by performing small affine transformations, keeping the 
volume and the number of particles in both phases con- 
stant. In the limit of an infinitesimal distortion of the 
system, Eq. (|2ip can be simplified as [36, ^] 



(22) 



where AE — Ei ^ Eq, as the probability distribution 
functions of the ensembles corresonding to the unper- 
turbed and the perturbed system converge, so that a 
separate sampling is no longer required. A third-order 
expansion in the inverse temperature (53l | 



AA 



{AE) (AE^) - {AEf 



kT 2(fcT)2 
(AE^) - 3 (AE^) {AE) 



2 {AE)'' 



6{kT)3 



(23) 



can be used to increase the precision of the simulation 
results [25. '3^ . The surface tension is then immediately 
obtained from AA/ Aa, since the distortion of the inter- 
face itself (as opposed to its increase in area) makes a 
negligible contribution to the free energy difference Il5| . 

In analogy with the Widom test-particle method |54| . 
this implementation of the variational route is also called 
the test-area method [H, HI] . Following Sampayo et al. 
(2^ . it can be applied to curved interfaces, where the 
afhne transformation scales one of the cartesian axes by 
the factor 1/(1-)-^) and the remaining ones by (1 -1-^)^/^. 
For ^ > 0, this creates an oblate shape and the area of 
the surface of tension is increased by |56| 



with the ellipticity of the average equimolar surface in 
the perturbed system given hy S = [1 — (1 -I- 
In the prolate case (^ < 0), the corresponding term is 
=■ [1 - (1-0"^]^^^ with |5i| 



Aa / arcsin r; 

^"H^(i-e)i/2 



o 



SAa 



(25) 



It can be shown that the first-order term in Eq. is 
equivalent to the Kirkwood-Buff |57| rnechanical route 
expression for the surface tension |58| . The higher- 
order terms therefore presumably capture the deviation 
between the mechanical and variational routes due to 
fluctuations or, equivalently, the contribution of non- 
equilibrium configurations to 7. Thus, the higher-order 
contribution to Eq. (|23p may be related to the closed ex- 
pression derived by Percus et al. (59| for the deviation 
between the actual free energy and an approximation 
based on the local pressure. 

From this point of view, the following aspects of the 
method merit further consideration: 

• While finite differences of higher order are taken 
into account for the energy, no such terms are con- 
sidered for the surface area here. Clearly, the vari- 
ance of AE is partly caused by the variance of Aa. 
The use of Rp for defining the surface area, cf. Eqs. 
and may lead to further deviations. 



• The variance of AE accounts for surface oscillations 
such as long wave-length capillary waves, which dir- 
ectly relate to equilibrium properties of the inter- 
face and therefore do not depend on the statistical 
mechanical ensemble ^|. However, it can also be 
influenced by fluctuations regarding p' (at constant 
V) or V (at constant p'). These modi are ensemble 
dependent, since they are coupled to the density of 
the vapour phase. Canonically, their amplitude in- 
creases with the total volume and is ill-defined in 
the thermodynamic limit V ^ 00. Therefore, the 
surface tension from the variational route may de- 
pend on the constraints imposed on the system by 
the ensemble. 

• Although the volume associated with each of the 
phases is invariant for test-area transformations, 
there is still a distortion of the sample with respect 
to the equilibrium conformation. The method is 
therefore limited to isotropic phases, since shearing 
an anisotropic phase will induce an elastic contri- 
bution to A^ from the bulk region as well. 



III. DEVIATION OF THE EQUIMOLAR 
RADIUS FROM CAPILLARITY 



From the Tolman equation in its approximate poly- 
nomial form, cf. Eq. pOp . the excess equimolar radius i] 



7 



can be related to the Tolman length 5 by 



= -S + 0{R-'), (26) 

so that its magnitude in the zero-curvature limit is ob- 
tained as 



Vo = -So 



(27) 



which is essentially equivalent to Eq. 

Both in the planar limit and in the presence of 
curvature effects, it is therefore possible to express the 
Tolman relations in terms of the easily accessible quant- 
ities 7] and if, rather than S and 1/R-y. The point of 
departure for such an expression is the exact closed form 
of the Tolman equation, cf. Eq. It should be re- 

called that this expression is derived from the Gibbs- 
Duhem equation, the Young-Laplace equation and the 
Gibbs adsorption equation [17| ; hence, it is based entirely 
on an axiomatic thermodynamic treatment. As opposed 
to truncated power series of the form of Eq. ([TU)) . the 
alternative description remains valid when the radius R^ 
becomes similar or smaller in magnitude than the Tol- 
man length. Polynomial expansions in terms of S / R-y 
necessarily fail to capture this limit. 
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Figure 2. Density profiles from canonical MD simulations of 
LJTS liquid drops at T = 0.75 e/k with equimolar radii of 

Rp = 9.977 ± 0.001 ( ), 12.029 ± 0.003 (- -), 13.974 ± 

0.002 o- (■ ■ ■ ) and 15.967 ± 0.001 a (— ), cf. Tabs. Hand HIl 



From the Young-Laplace equation, it follows that 



dR^ 

dip 



1 d7 

ip dip 



J_ 

ip^ 



(28) 



while the reduced length scale appearing in the Tolman 
equation can be transformed to 



5 

R-y 



+ 70 
7 



1, 



(29) 



by using Eqs. (H]), (g]), ([5]) and HT])- The Tolman relation 
can thus be converted to 



(30) 



c?7 


_2y 




■ 5 ' 


' 1 


■ 5 ' 




dip 




R^ ^ 


R~f 


+ 3 


R^ 






+ 70 



7 



(31) 



This representation of the Tolman result is fully equival- 
ent to Eq. (O. 

For ip> 0, further considerations are required. There, 
the curvature dependence of 7 as specified by Eq. ((3T|) 
is only self-consistent under an additional condition. To 
demonstrate this, it is helpful to consider the exact Tol- 
man equation in a different form 



^ Q ['^ ^ '^^'^^] ^ 7o??[7o + Vf] 



(32) 



which follows from Eq. pip by expanding the cubic term 
where ^ has been defined as 



c = 



7 



70^ 



(33) 



For the sake of brevity, the notation qt = lini;^_i.o d''q/dip^ 
is used here for the z-th derivative of a quantity q in the 
zero-curvature limit. The slope of 7 can be obtained by 
inserting 



37o7i, 



Co = (7^), 
into Eq. ([5^ . which yields 

71 = 2ryo- 

Expanding the excess equimolar radius as 

V ^vo + mv + oiip'^), 



(34) 



(35) 



(36) 



and inserting this expression as well as Eq. (I35p into the 
planar limit for Eq. (15^ leads to 



(7^)2 = 1270%', 



and 



7o72 



(37) 



(38) 



It is by considering the zero-curvature limit for the third 
derivative of 7"^ that a theorem for the slope of rj can 
now be deduced. Based on Eqs. and (P7)) . a Taylor 
expansion for d{'j^)/dip in terms of ip 

+ + Oiip% (39) 



yields 



V Jq \dip 
670*70 + 670770 <y5 4 



0{ip3). (40) 
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From Eqs. to (HP)) 



670 {lorn + rfo) + ^ (7')3 = + O(^), (41) 



follows by applying the full Tolman equation, cf. Eq. 
in the planar limit. However, this implies 

m = — , 

7o 



(42) 



which constitutes a necessary boundary condition for the 
Tolman approach in terms of 77 and ip. 

Thus, while there is a direct correspondence between 
i5o and 770, no such relation exists in case of Si and 771, 
i.e., the respective derivatives (in terms of ip) in the zero- 
curvature limit; instead, 771 is fully determined by 770 and 
thus by So, the Tolman length of the planar interface. 
This means that data on the excess equimolar radius for 
large radii have a double significance regarding the planar 
limit: on the one hand, they can be extrapolated to f — 
0, leading to an estimate for the planar Tolman length 
and the curvature dependence of 7 to first order in terms 
of (f or 1/ Rj] on the other hand, the slope of rj is in itself 
relevant, since its zero-curvature limit rji also provides 
information on rjQ. 

The equivalent of the exact Tolman equation in terms 
of the excess equimolar radius 77 and the pressure differ- 
ence characterized by ip is Eq. (|3ip . An expansion as a 
power series, analogous to Eq. PU)) . can be expressed as 

7 = 7o + 277o(^-^(^2 + 0((^3). (43) 
70 

The planar limit, where higher order terms can be neg- 
lected, can be treated accurately with expressions like 
Eq. ([43|) . Away from the planar limit, Eq. ([3T|) applies 
without any further condition (since the boundary con- 
dition for the slope of 77 is only relevant for (p ^ 0), while 
Eq. ([43]) becomes an approximation. 



IV. 



THE EXCESS EQUIMOLAR RADIUS FROM 
MOLECULAR SIMULATION 



With the mardyn MD progra m, developed by 
Bernreuther and co-workers [GOUS^I, the canonical en- 
semble was simulated for small systems, corresponding 
to equilibrium conditions for nanoscopic liquid drops 
surrounded by supersaturated vapours. The truncated- 
shifted Lennard- Jones (LJTS) pair potential 




Ushift, 



for r < Vc 
for r > Tc 



(44) 



with the size parameter a, the energy parameter e and 
a cutoff at Tc = 2.5 cr is applied as a fluid model here, 
including a shift by u^Mit = 4e [{(j/rc)^ - (cr/r-c)^^] to 
make the potential continuous. The LJTS model is an ad- 
equate basis for investigating bulk and interfacial proper- 
ties of simple spherical conformal fluids (e.g., noble gases 



and methane) on a molecular level, cf. Vrabec et al. |25| . 
On account of this, numerous studies on nanoscopic li- 
quid drops have been reported [11 [M HI I63l467j. The 
LJTS fluid can thus be regarded as a key benchmark for 
theoretical and simulation approaches to the problem of 
curved vapour-liquid interfaces. 

Certain of the general properties of this simple model, 
taking only short-range interactions into account, can be 
assumed to carry over to polar fluids as well (68j , except 
for temperatures in the vicinity of the critical point. It 
is clear, however, that a qualitatively different behaviour 
should be expected for liquid drops formed by water with 
and without ionic species (69l . [70j , liquid crystals fn\ and 
similar complex organic molecules. Such systems are bey- 
ond the scope of the present study. 

Liquid drops are investigated at temperatures between 
T = 0.65 and 0.95 e/k, covering most of the range 
between the triple point temperature (which is « 0.55 
according to Bolhuis and Chandler |72i], w 0.618 as de- 
termined by Toxvaerd [t^ and ~ 0.65 e/k according to 
van Meel et al. [65|) and the critical temperature which 
several independent studies have consistently obtained 
as 1.08 e/k for the LJTS fluid [H III, The Ver- 

let leapfrog algorithm is employed to solve the classical 
equations of motion numerically with an integration time 
step of 0.002 in Lennard- Jones time units, i.e., ay^m/e, 
where m is the mass of a particle. Cubic simulation 
volumes with 290 to 126 000 particles, applying the peri- 
odic boundary condition, are equilibrated for at least 2 
000 time units. Subsequently, spherically averaged dens- 
ity profiles p{z), with their origin (z = 0) at the centre 
of mass of the whole system, are constructed with a bin- 
ning scheme based on equal volume concentrical spheres 
using sampling intervals between 1 000 and 40 000 time 
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Figure 3. Density profiles from canonical MD simulations of 
LJTS liquid drops at T = 0.65 and 0.95 e/k, showing the av- 
erage densities from simulation (•) and exponential approx- 

imants ( ). The steeper profile corresponds to the lower 

temperature. 
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Figure 4. Density profiles from a single canonical MD sim- 
ulation of a LJTS liquid drop at T = 0.85 e/k, showing the 
average densities from simulation (o) and exponential approx- 
imants (lines) corresponding to the sampling intervals 2 000 - 

3 000 (• • ■ ; green), 3 000 - 4 000 (• • • ; red), 4 000 - 5 000 

( ; blue) and 5 000 - 6 000 time units ( — ; black) after the 

onset of the simulation. The standard deviation between the 
densities at an infinite distance from the interface, according 
to the exponential fits for all sampling intervals of a single 
MD simulation, is used to determine the error of the bulk 
densities here. 
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Figure 5. Equimolar radius Rp as a function of the capillarity 
radius Rn for LJTS liquid drops, from density profiles and 
bulk pressures determined with canonical MD simulations at 
T = 0.75 (□) and 0.85 e/k (o), in comparison with results 
from previous work of Vrabec et al. [25| at T = 0.75 (■) and 
0.85 e/k (•), using pressure differences based on evaluating 
the IK tensor in the (approximately) homogeneous regions 
inside and outside the liquid drop. The continuous diagonal 
line is defined by Rp — i?„ and thus corresponds to an excess 
equimolar radius of 77 = 0, while the dotted lines correspond 
to = ±0.5 a. 



units, depending on the (expected) total simulation time, 
to gather multiple samples for each system. Examples of 
the density profiles obtained according to this method 
are shown in Figs. [2]- 21 

The density profiles of LJTS vapour-liquid interfaces 
are known to agree well with an expression based on two 
hyperbolic tangent terms, to which p{z) has been suc- 
cessfully correlated for liquid drops by Vrabec et al. . 
The present method, however, merely requires the bulk 
densities p' and p" corresponding to a certain value of p 
or which are determined here by correlating the outer 
part of the density profile and extrapolating it to regions 
far from the interface. The densities of the coexisting 
fluid phases are thus deduced from simulation results by 
correlating the exponential terms 

p' = p{z) + a' exp {I3'[z - z']) , 
p" = p(z)-a"exp(/3"[^"-z]), (45) 

to the data for the inner- and outermost spherical bins 
of the density profiles, cf. Fig. [3l These terms, which are 
based on those employed by Lekner and Henderson [5S^ . 
asymptotically agree with the hyperbolic tangent expres- 
sion of Vrabec et al. (25l |. From the liquid and vapour 
densities p' and p" of the fit to Eq. (j45|) . the equimolar 
radius Rp is calculated according to Eq. ^ . The respect- 
ive margins of error are obtained as standard deviations 
from the profiles belonging to different sampling intervals 
of the same MD simulation, cf. Fig. 01 of which there are 
at least three in all cases. The corresponding pressures 



p' and p" are computed by canonical MD simulation of 
the bulk fluid at the respective densities. 

For the surface tension in the zero-curvature limit, 
the values 7o(0.65e/fc) = 0.680 ± 0.009, 7o(0.75£//c) = 
0.493 ± 0.008, 7o(0.85e/A:) = 0.317±0.007 and 7o(0.95 e/ 
k) = 0.158 ± 0.006 ecr"^ are taken from the correlation 
of Vrabec et al. [1^ ; the error corresponds to the indi- 
vidual data points for 70 from the same source. In case 
of T = 0.9 e/k, the higher precision of the computations 
of van Giessen and Blokhuis [24] is exploited, using the 
value 7o = 0.227 ± 0.002 ecr"^ obtained from a linear 
fit to data for the curved interface [l^], cf. Fig. [1] The 
assumption made for the error is rather generous in this 
case, considering the even higher confidence suggested 
by the agreement between the individual data points for 
(pRp. 

Combining these quantities leads to the capillarity ra- 
dius i?K and the excess equimolar radius 77. Note that the 
margin of error for 77, as indicated in Tab.U contains con- 
tributions quantifying the accuracy of 70 and the preci- 
sion of the MD simulations of the liquid drop itself as well 
as those of the homogeneous vapour and liquid phases. 
While the vapour pressure p" and the equimolar radius 
Rp could be obtained with a high precision, the liquid 
pressure and the surface tension in the zero-curvature 
limit are major sources of uncertainty here. In both cases, 
methodical changes can be expected to increase the pre- 
cision significantly: regarding 70, it can be seen from Fig. 
[Hthat it is now possible to reach a level of confidence bey- 
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Table I. An analysis of the error of the excess equimolar radius rj of LJTS liquid drops at the temperature T = 0.75 e/k. 
The number of particles A'', the volume V of the periodic simulation box and the total simulation time t for the simulations 
of the liquid drops are indicated alongside the contributions to the uncertainty of rj from the pressure p' of the liquid phase 
(determined by canonical MD simulation of the bulk liquid), the surface tension 70 of the planar vapour- liquid interface, cf. 
Vrabec et al. [25l |. the vapour pressure p" (analogous to p') and the equimolar radius Rp (from the density profiles of the liquid 
drops). Note that the time unit, i.e., a^/rnje, corresponds to 500 simulation time steps here. The flat symbols (b) indicate the 
fraction of the margin of error for rj due to the respective quantities. All values are given in Lennard- Jones units, and the error 
in terms of the last digit is specified in parentheses. In the subsequent discussion, the cases where the uncertainty of rj exceeds 
a are disregarded. 
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0.16(1) 


82% 
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0.5% 
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6 


000 


0.14(3) 
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0.113(2) 
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0.01042(4) 
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0.28% 


0.3(3) 


5 161 
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64 
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0.096(3) 


63% 


34% 


0.0104(1) 


2.4% 
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0.82 % 
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0.090(2) 
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40% 
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0.75% 


12.029(3) 


0.57% 
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10 241 
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0.58 % 
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0.29 % 
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32 
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66% 
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0.00974(4) 


1.2% 


14.981(6) 


0.78% 
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15 237 
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0.35 % 
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17 113 


169 418 


6 
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16.689(4) 


0.18% 


2(2) 
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6 
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3(5) 



Table II. Number of particles A'^, volume V of the periodic simulation box and temperature T of the present canonical ensemble 
MD simulations of the LJTS fluid and equilibrium properties of the liquid drop as well as the surrounding vapour, i.e., the 
respective densities p', p" and pressures p' , p" as well as the capillarity radius i?^, the equimolar radius Rp and the excess 
equimolar radius rj. For drop radii above 8 a, these values can be reliably regarded as identical with those corresponding 
to the present theoretical approach, which is highlighted with the bold typeface. In case of smaller radii (cursive typeface), 
inaccuracies can arise due to the application of exponential approximants, cf. Fig. [3] and Eq. (|45ll . so that the respective values 
can, at present, be acknowledged as phenomenological quantities only. All values are given in Lennard- Jones units, and the 
error in terms of the last digit is specified in parentheses. 



iV 




T [e/k] 




P" 


P' [so-''] 


P" [ea-^[ 


R. \o] 


Rp [a] 


V H 


291 


8 999 


0.65 


0.857{5) 


0.0090(2) 


0.65(8) 


0.0054(1) 


2.1(3) 


3.90(1) 


1.8(3) 


1 022 


17 576 


0.65 


0.830(1) 


0.00651(7) 


0.22(2) 


0.00397(4) 


6.3(6) 


6.407(2) 


0.1(6) 


497 


10 648 


0.75 


0.81(1) 


0.0214(6) 


0.6(1) 


0.0135(4) 


1.8(5) 


4.33(5) 


2.5(5) 


1 418 


21 952 


0.75 


0.777(1) 


0.0173(1) 


0.16(1) 


0.01136(5) 


6.5(6) 


6.883(3) 


0.4(6) 


3 762 


39 304 


0.75 


0.7721(2) 


0.01566(6) 


0.113(2) 


0.01042(4) 


9.7(4) 


9.977(1) 


0.3(4) 


5 161 


54 872 


0.75 


0.7703(2) 


0.0156(2) 


0.096(3) 


0.0104(1) 


11.5(5) 


11.089(4) 


-0.5(6) 


6 619 


74 088 


0.75 


0.7697(2) 


0.01506(4) 


0.091(2) 


0.01007(2) 


12.3(5) 


12.029(3) 


-0.2(5) 


10 241 


110 592 
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0.075(2) 


0.00974(4) 


15.2(8) 


14.981(6) 


-0.2(8) 


15 237 


166 375 
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0.662(2) 


0.0825(2) 


0.169(7) 
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ond that of the data of Vrabec et al. [l^ which are also 
used here. For the pressure of the hquid p' , approaches 
based on the chemical potential, which can be determ- 
ined in any region of the simulation volume (including 
the vapour phase) , can be expected to lead to significant 
improvements in combination with a reliable equation of 
state or high-precision simulations in the grand canonical 
ensemble. 

A full summary of the simulation results where could 
be determined with error bars smaller than a is given in 
Tab. ini Note that to achieve full consistency with the 
Tolman approach, the bulk densities p' and p" from Eq. 
(|45| have to match those of the bulk fluid at the same 
temperature and chemical potential as the two-phase sys- 
tem. Regarding liquid drops with Rp > 8 a, this is cer- 
tainly the case, since constant density regions coexisting 
with the interface are actually present, cf. Fig. How- 
ever, the values determined for the smallest drops here 
rely on the validity of the correlation given by Eq. ((45|) 
and can be considered valid only as far as this expression 
itself does not introduce any major deviations, an asser- 
tion that remains open to further examination; a version 
of the present method computing p' via p could resolve 
this issue. 



V. DISCUSSION 

Previous authors have made qualitatively contradict- 
ory claims on the magnitude of the Tolman length as 
well as its sign: Tolman himself expected 5 to be pos- 
itive and smaller than the length scale of the dis per sive 
interaction, a conjecture that Kirkwood and Buff [57j af- 
firmed from a statistical mechanical point of view, based 
on a mechanical approach. Subsequent studies, however, 
have also found S to be negligible or even equal to zero 
(20I m, [7^, positive and larger than a [l^ls^, negat- 
ive with -a < S < [22!, [3, or negative and di- 
verging {60 — —00) in the planar limit [78|, while others 
have claimed that the sign of S is curvature dependent it- 
self [7^ [s^l ■ Thereby, they have only proven the mutual 
inconsistency of their assumptions and methods, while 
nothing is truely known about S and the dependence of 
the surface tension on curvature. 

The new approach introduced in Section IIIII is strictly 
based on axiomatic thermodynamics and relies on the 
fact that i5o = —770 holds in the planar limit. From the 
values for rj reported in bold face in Tab. [TTl correspond- 
ing to Rp > 8 (J, the excess equimolar radius for liquid 
drops of the LJTS fluid is unequivocally shown to be 
smaller in magnitude than cr/2, while its remains unclear 
whether it is positive, negative, of both signs (depend- 
ing on the curvature) or equal to zero. Since this means 
that at the present level of accuracy, no significant de- 
pendence of 7 on the radius of the liquid drop could be 
detected, the statement of Mareschal et al. [8l| regarding 
cylindrical interfaces also applies here: considering «the 
large fluctuations in the bulk liquid phase» , cf. the error 



analysis presented in Tab. HI «we tentatively conclude 
that the surface tension is independent of the curvature 
of the liquid-vapor interface or else that this dependence 
is very weak.» 

The only view that can be deflnitely dismissed is that 
of a large and positive Tolman length, previously held by 
some of the present authors on the basis of results from 
the mechanical route to the surface tension, employing 
the IK pressure tensor [2^|66|. As Fig. [S] shows, the pre- 
vious simulation results are actually consistent with those 
from the present study if they are interpreted in terms 
of the radii R^ and Rp. Thereby, following the approach 
of van Giessen and Blokhuis ^] , only the density profile 
and the pressure in the homogeneous regions inside and 
outside the liquid drop are taken into account, whereas 
the normal pressure along the interface is not considered 
at all. Since the deviation between present and previous 
data disappears in such a representation, the disagree- 
ment must be caused by the inadequacy of the pressure- 
tensor (mechanical) route implemented by Thompson et 
al. [Ill, as pointed out by Henderson [39|, Hlj. Possible 
sources of error for this approach are outlined in Section 
im Nonetheless, more detailed methodological invest- 
igations are expedient to determine which approxima- 
tions are actually responsible for major inaccuracies, and 
whether they can be corrected or whether the pressure- 
tensor route to the surface tension has to be discarded 
altogether. 
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